1 Introduction

The objectives of this notebook are:

  1. Include the scores computed in the previous notebook to the metadata of the Seurat object.
  2. Project such scores and annotations in the UMAP to get a sense of where our predicted doublets are located.
  3. Exclude the doublets and save the object.

2 Pre-processing

2.1 Load packages

library(AnnotationDbi)
library(org.Hs.eg.db)
library(Seurat)
library(Matrix)
library(tidyverse)

2.2 Parameters

# Paths
path_to_obj <- here::here("scRNA-seq/results/R_objects/seurat_integrated_with_doublets.rds")
path_to_doubl_annot <- here::here("scRNA-seq/2-QC/5-batch_effect_correction/1-doublet_exclusion/tmp/doublet_final_annotations_and_pDNN.rds")
path_to_save <- here::here("scRNA-seq/results/R_objects/seurat_without_doublets_missing_reintegration.rds")


# Functions
source(here::here("scRNA-seq/bin/utils.R"))

2.3 Load data

tonsil <- readRDS(path_to_obj)
doublet_annot <- readRDS(path_to_doubl_annot)

3 Add metadata

if (all(rownames(doublet_annot) == colnames(tonsil))) {
  tonsil$pDNN_hashing <- doublet_annot$pDNN_hashing
  tonsil$pDNN_scrublet <- doublet_annot$pDNN_scrublet
  tonsil$pDNN_union <- doublet_annot$pDNN_union
} else {
  warning("barcodes are not equal")
}

4 Visualize

4.1 Proportion of doublet nearest neighbors (pDNN)

pDNN_vars <- c("pDNN_hashing", "pDNN_scrublet", "pDNN_union")
pDNN_ggs <- purrr::map(pDNN_vars, function(x) {
  p <- feature_plot_doublets(seurat_obj = tonsil, feature = x)
  p
})
names(pDNN_ggs) <- pDNN_vars
pDNN_ggs
## $pDNN_hashing

## 
## $pDNN_scrublet

## 
## $pDNN_union

4.2 Doublet annotations

doublet_vars <- c("HTO_classification.global", "scrublet_predicted_doublet",
                  "has_high_lib_size")
doublet_ggs <- purrr::map(doublet_vars, function(x) {
  p <- DimPlot(tonsil, group.by = x, pt.size = 0.3) +
    ggtitle(x) +
    theme(
      plot.title = element_text(size = 14, hjust = 0.5),
      axis.title = element_text(size = 13),
      axis.text = element_text(size = 11),
      legend.text = element_text(size = 11)
    )
  p
})
names(doublet_ggs) <- doublet_vars
doublet_ggs
## $HTO_classification.global

## 
## $scrublet_predicted_doublet

## 
## $has_high_lib_size

# Zoom in for cell hashing
umap_hashing <- tonsil@reductions$umap@cell.embeddings %>%
  as.data.frame() %>%
  mutate(HTO_classification.global = tonsil$HTO_classification.global) %>%
  ggplot(aes(UMAP_1, UMAP_2, color = HTO_classification.global)) +
    geom_point(size = 0.01, alpha = 0.5) +
    facet_wrap(~HTO_classification.global) +
    theme_classic()
umap_hashing

4.3 Canonical markers

canonical_markers <- c("CD79A", "CD79B", "CD3D", "CD3E", "NKG7", "LYZ", "FDCSP")
canonical_markers_gg <- purrr::map(canonical_markers, function(x) {
  p <- feature_plot_doublets(seurat_obj = tonsil, feature = x)
  p
})
names(canonical_markers_gg) <- canonical_markers
canonical_markers_gg
## $CD79A

## 
## $CD79B

## 
## $CD3D

## 
## $CD3E

## 
## $NKG7

## 
## $LYZ

## 
## $FDCSP

4.4 pre-Plasmablast and tingible-body macrophages

We want to be particularly careful with two cell types: pre-plasmablasts and tingible-body macrophages, since they might be confused by doublets (specially using scrublet). Thus, we will visualize a prePB and apoptosis signature:

# Define signatures
apoptosis_signature <- AnnotationDbi::select(
  org.Hs.eg.db,
  keys = "GO:0006915",
  keytype = "GO",
  columns = "SYMBOL"
)$SYMBOL
apoptosis_signature <- unique(apoptosis_signature)
prePB_signature <- c("RASSF6", "FRZB", "HOPX", "BRNL9", "FGFR1")
signatures_list <- list(
  prePB_signature = prePB_signature,
  apoptosis_signature = apoptosis_signature
)


# Compute score
tonsil <- AddModuleScore(
  tonsil,
  features = signatures_list,
  name = c("prePB_signature", "apoptosis_signature")
)


# Visualize
prePB_gg <- feature_plot_doublets(
  seurat_obj = tonsil,
  feature = "prePB_signature1"
)
apopotosis_gg <- feature_plot_doublets(
  seurat_obj = tonsil,
  feature = "apoptosis_signature2"
)
prePB_gg

apopotosis_gg

5 Exclude doublets

tonsil
## An object of class Seurat 
## 28504 features across 347262 samples within 1 assay 
## Active assay: RNA (28504 features, 3000 variable features)
##  3 dimensional reductions calculated: pca, umap, harmony
table(tonsil$HTO_classification.global)
## 
## Doublet      NA Singlet 
##   80577   50072  216613
tonsil <- subset(tonsil, subset = HTO_classification.global != "Doublet")
tonsil
## An object of class Seurat 
## 28504 features across 266685 samples within 1 assay 
## Active assay: RNA (28504 features, 3000 variable features)
##  3 dimensional reductions calculated: pca, umap, harmony

6 Save

saveRDS(tonsil, path_to_save)

7 Session Information

sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
## 
## Matrix products: default
## BLAS:   /apps/R/3.6.0/lib64/R/lib/libRblas.so
## LAPACK: /home/devel/rmassoni/anaconda3/lib/libmkl_rt.so
## 
## locale:
##  [1] LC_CTYPE=C                 LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.5.0        stringr_1.4.0        dplyr_1.0.4          purrr_0.3.4          readr_1.3.1          tidyr_1.1.0          tibble_3.0.1         ggplot2_3.3.0        tidyverse_1.3.0      Matrix_1.2-18        Seurat_3.2.0         org.Hs.eg.db_3.10.0  AnnotationDbi_1.48.0 IRanges_2.19.14      S4Vectors_0.24.4     Biobase_2.45.1       BiocGenerics_0.32.0  BiocStyle_2.14.4    
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1          backports_1.1.7       plyr_1.8.6            igraph_1.2.5          lazyeval_0.2.2        splines_3.6.0         listenv_0.8.0         digest_0.6.20         htmltools_0.4.0       viridis_0.5.1         fansi_0.4.1           magrittr_1.5          memoise_1.1.0         tensor_1.5            cluster_2.1.0         ROCR_1.0-11           globals_0.12.5        modelr_0.1.8          colorspace_1.4-1      blob_1.2.1            rvest_0.3.5           rappdirs_0.3.1        ggrepel_0.8.2         haven_2.3.1           xfun_0.14             crayon_1.3.4          jsonlite_1.7.2        spatstat_1.64-1       spatstat.data_1.4-3   survival_3.1-12       zoo_1.8-8             ape_5.3               glue_1.4.1            polyclip_1.10-0       gtable_0.3.0          leiden_0.3.3          future.apply_1.5.0    abind_1.4-5           scales_1.1.1          DBI_1.1.0             miniUI_0.1.1.1        Rcpp_1.0.6            viridisLite_0.3.0     xtable_1.8-4          reticulate_1.16       bit_1.1-15.2          rsvd_1.0.3            htmlwidgets_1.5.1     httr_1.4.2            RColorBrewer_1.1-2    ellipsis_0.3.1        ica_1.0-2             farver_2.0.3          pkgconfig_2.0.3      
##  [55] uwot_0.1.8            dbplyr_1.4.4          deldir_0.1-25         here_0.1              labeling_0.3          tidyselect_1.1.0      rlang_0.4.10          reshape2_1.4.4        later_1.0.0           munsell_0.5.0         cellranger_1.1.0      tools_3.6.0           cli_2.0.2             generics_0.0.2        RSQLite_2.1.2         broom_0.5.6           ggridges_0.5.2        evaluate_0.14         fastmap_1.0.1         yaml_2.2.1            goftest_1.2-2         knitr_1.28            bit64_0.9-7           fs_1.4.1              fitdistrplus_1.1-1    RANN_2.6.1            pbapply_1.4-2         future_1.17.0         nlme_3.1-148          mime_0.9              xml2_1.3.2            compiler_3.6.0        rstudioapi_0.11       plotly_4.9.2.1        png_0.1-7             spatstat.utils_1.17-0 reprex_0.3.0          stringi_1.4.6         lattice_0.20-41       vctrs_0.3.6           pillar_1.4.4          lifecycle_0.2.0       BiocManager_1.30.10   lmtest_0.9-37         RcppAnnoy_0.0.16      data.table_1.12.8     cowplot_1.0.0         irlba_2.3.3           httpuv_1.5.3.1        patchwork_1.0.0       R6_2.4.1              bookdown_0.19         promises_1.1.0        KernSmooth_2.23-17   
## [109] gridExtra_2.3         codetools_0.2-16      MASS_7.3-51.6         assertthat_0.2.1      rprojroot_1.3-2       withr_2.4.1           sctransform_0.2.1     mgcv_1.8-31           hms_0.5.3             grid_3.6.0            rpart_4.1-15          rmarkdown_2.2         Rtsne_0.15            shiny_1.4.0.2         lubridate_1.7.8